In [12]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 5, 101)
s = np.sin(t**2) * np.exp(-(t-t[50])**2)

lm = ['.', '-', '--', ':', '-.']

plt.figure(figsize=(12, 5))
for i in range(len(lm)):
    plt.plot(t, s+i, lm[i], linewidth=5, label='line style: ' + lm[i])
plt.legend(fontsize=25)
plt.xlim(-0.25, 8.5)
plt.show()
In [56]:
import seaborn as sns
import matplotlib.pyplot as plt
 
x = np.arange(0., np.e, 0.01)
y1 = np.exp(-x)
y2 = np.log(x)
 
fig = plt.figure(figsize=(8, 6))
ax1 = fig.add_subplot(111)
ax1.plot(x, y1)
ax1.set_ylabel('Y values for exp(-x)', fontsize=20)
ax1.set_title("Double Y axis", fontsize=20)
 
ax2 = ax1.twinx()  # this is the important function
ax2.plot(x, y2, 'r')
ax2.set_xlim([0, np.e])
ax2.set_ylabel('Y values for ln(x)', fontsize=20)
ax2.set_xlabel('Same X for both exp(-x) and ln(x)', fontsize=20)

plt.show()
D:\Anaconda\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in log
  
In [85]:
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(100, 100000, 1001)
y = np.log(x)
plt.figure(figsize=(10,7))
plt.subplot(311)
plt.semilogx(x, y, 'r', label='log(x) linear y')
plt.legend(fontsize=20)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.subplot(312)
plt.semilogy(x, y, 'g', label='log(y) linear x')
plt.legend(fontsize=20)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.subplot(313)
plt.loglog(x, y, 'b', label='log(x) log(y)')
plt.suptitle('Log axis', fontsize=25)
plt.legend(fontsize=20)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.show()
In [102]:
import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(10,8))
fig, ax = plt.subplots()
ax.plot(np.linspace(0,2*np.pi, 101), np.sin(np.linspace(0, 2*np.pi, 101)), 'r', linewidth=2)
for line in ax.get_xticklines() + ax.get_yticklines():
    line.set_markeredgewidth(10)
    line.set_markersize(5)
    line.set_markeredgecolor('b')
    plt.xticks(fontsize=18)
    plt.yticks(fontsize=18)
plt.suptitle('Special ticks', fontsize=25)
plt.show()
<Figure size 720x576 with 0 Axes>
In [83]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s = abs(np.sin(t))

plt.figure(figsize=(20, 5))
plt.subplot(121)
ax = plt.semilogy(t, s)
# Major grid lines.
plt.grid(b=True, which='major', color='r', alpha=0.5, linewidth=1,  visible=True, ls='--')
# Minor grid lines.
plt.grid(b=True, which='minor', color='g', alpha=0.5, linewidth=1,  visible=True, ls=':')
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)


plt.subplot(122)
plt.plot(t, s)
ax = plt.gca()
# Grid lines on Y axis.
ax.yaxis.grid(ls=':', color='r')
ax.set_yticks(np.linspace(-1, 1, 11))
# Grid lines on X axis.
ax.xaxis.grid(ls='--', color='g')
ax.set_xticks(np.linspace(0, 10, 11))
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.suptitle('Logarithm ticks and grid lines', fontsize=22)

plt.show()
In [67]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 101)
n = np.random.random(101)
plt.figure(figsize=(15, 6))
plt.subplot(121)
plt.stem(t, n, markerfmt='*', linefmt=':')

plt.subplot(122)
t = np.linspace(0, 10, 101)
n = np.random.random(101)
plt.stem(t, n, markerfmt='^', linefmt='-')
plt.suptitle('stem', fontsize=28)
plt.show()
In [41]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 5, 51)
s = np.sin(t**2) * np.exp(-(t-t[25])**2)
markers = ['.', 'o', 's', 'P', 'X', '*', 'p', 'D', '<', '>', '^',
           'v', '1', '2', '3', '4', '+', 'x', '|', '_']
colors = ['r', 'g', 'b', 'c', 'k', 'm', 'y']
n = len(markers)
l = n // 3
if l*3 < n:
    l += 1

plt.figure(figsize=(15, 25))
for i in range(n):
    index1 = int(np.random.random()*7)
    index2 = int(np.random.random()*7)
    plt.subplot(l, 3, i+1)
    plt.plot(t, s, color='#555555', linewidth=1)
    plt.scatter(t, s, marker=markers[i], s=45, facecolor=colors[index1], edgecolor=colors[index2], linewidth=1, label=markers[i])
    plt.legend(fontsize=15)

plt.subplot(l, 3, i+2)
plt.scatter(t, s, marker=7, s=25, facecolor=colors[index1], edgecolor=colors[index2], linewidth=1)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.xlabel('Time [s]', fontsize=20)
#plt.ylabel('Amplitude', fontsize=20)
plt.legend(fontsize=15)
plt.suptitle('Marker styles', fontsize=25)
plt.show()
No handles with labels found to put in legend.
No handles with labels found to put in legend.
In [14]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(-3, 3, 601)
X, Y = np.meshgrid(t, t)
Z = np.sin(X**2+Y**2) * np.exp(-X**2-Y**2)

cmaps_str='Accent, Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, \
CMRmap, CMRmap_r, Dark2, Dark2_r, GnBu, GnBu_r, Greens, Greens_r, Greys, Greys_r, \
OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, Pastel1, Pastel1_r, \
Pastel2, Pastel2_r, PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, \
PuRd_r, Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r, RdPu, RdPu_r, RdYlBu, RdYlBu_r, \
RdYlGn, RdYlGn_r, Reds, Reds_r, Set1, Set1_r, Set2, Set2_r, Set3, Set3_r, Spectral, Spectral_r, \
Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r, YlOrRd, YlOrRd_r, afmhot, \
afmhot_r, autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cividis, \
cividis_r, cool, cool_r, coolwarm, coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, flag, \
flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, \
gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, \
gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, hot, hot_r, hsv, hsv_r, inferno, inferno_r, \
jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, pink, pink_r, plasma, \
plasma_r, prism, prism_r, rainbow, rainbow_r, seismic, seismic_r, spring, spring_r, summer, summer_r, \
tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, viridis, viridis_r, winter, winter_r'

cmaps = cmaps_str.split(',')
n = len(cmaps)
m = n // 10
if m*10 < n:
    m += 1    
plt.figure(figsize=(25, 50))
for i in range(n):
    plt.subplot(m, 10, i+1)
    cmap = cmaps[i].strip()
    plt.pcolormesh(X, Y, Z, cmap=cmap)
    plt.title(cmap)
    
plt.show()
In [60]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

t = np.linspace(-3, 3, 601)
X, Y = np.meshgrid(t, t)
Z = np.sin(X**2+Y**2) * np.exp(-X**2-Y**2)

cdict = {'red': ((0., 1, 1),
                 (0.05, 1, 1),
                 (0.11, 0, 0),
                 (0.66, 1, 1),
                 (0.89, 1, 1),
                 (1, 0.5, 0.5)),
         'green': ((0., 1, 1),
                   (0.05, 1, 1),
                   (0.11, 0, 0),
                   (0.375, 1, 1),
                   (0.64, 1, 1),
                   (0.91, 0, 0),
                   (1, 0, 0)),
         'blue': ((0., 1, 1),
                  (0.05, 1, 1),
                  (0.11, 1, 1),
                  (0.34, 1, 1),
                  (0.65, 0, 0),
                  (1, 0, 0))}

my_cmap = matplotlib.colors.LinearSegmentedColormap('map',cdict,256)

plt.figure(figsize=(15, 4.8))
plt.subplot(131)
plt.pcolormesh(X, Y, Z, cmap=my_cmap)
plt.title('Create your colormap', fontsize=15)
plt.subplot(132)
plt.pcolormesh(X, Y, Z, cmap='jet')
plt.subplot(133)
plt.pcolormesh(X, Y, Z, cmap='RdGy_r')
plt.show()
In [73]:
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

x = np.linspace(-10, 10, 2001)
y = x
[X, Y] = np.meshgrid( x, y )
Z = np.exp( -(X**2+Y**2) )

# 自定义colormap.
cm = mpl.colors.LinearSegmentedColormap.from_list('cmap', ['#333333', '#660588', '#98F5FF', '#00FF00', '#FFFF00','#FF0000', '#8B0000'], 256)
plt.figure(figsize=(5, 5))
plt.pcolormesh( X, Y, Z, cmap=cm )
plt.title('Make youe own color map', fontsize=23)
plt.show()
In [90]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal._savitzky_golay import savgol_filter as sgolay
from wavelets import WaveletAnalysis as WA
from obspy import read

st = read()
tr = st[0]
dt = tr.stats.delta
print('dt = ', dt)
fs = 1 / dt

tr.detrend()
tr.data = tr.data - np.mean(tr.data)
tr.taper(0.025)
f1 = 1
f2 = 10
tr.filter('bandpass', freqmin=f1, freqmax=f2, corners=4, zerophase=True)
t = np.arange(tr.stats.npts) * dt
s = tr.data

s = s /max(abs(s))

ss = sgolay(s, n//60*2+1, 3)
#ss = s

wa = WA(s, dt=dt)
p = (wa.wavelet_power)
tt = wa.time
ff = 1 / wa.scales

re_s = wa.reconstruction()

plt.figure(figsize=(18, 15))
plt.subplot(311)
plt.plot(t, ss, linewidth=1, label='smoothed signal')
plt.plot(t, s, linewidth=4, label='signal')
plt.plot(t, re_s, linewidth=1, label='reconstructed signal')


plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.legend(fontsize=20)
plt.xlim(t[0], t[-1])

plt.subplot(312)
plt.pcolormesh(tt, ff, p, cmap='jet')
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Frequency [Hz]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.ylim(f1, f2)

ax = plt.subplot(313)
tr.spectrogram(log=False, axes=ax, cmap='jet')
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Frequency [Hz]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.ylim(f1, f2)
plt.ylim(f1, f2)
plt.show()
dt =  0.01
D:\Anaconda\lib\site-packages\obspy\signal\detrend.py:31: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if not np.issubdtype(data.dtype, float):
D:\Anaconda\lib\site-packages\obspy\core\trace.py:2111: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  if not np.issubdtype(self.data.dtype, float):
D:\Anaconda\lib\site-packages\scipy\signal\_arraytools.py:45: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  b = a[a_slice]
D:\Anaconda\lib\site-packages\wavelets\transform.py:104: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  wavelet_data[slices],
D:\Anaconda\lib\site-packages\scipy\fftpack\basic.py:160: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  z[index] = x
D:\Anaconda\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
In [17]:
import numpy as np
import matplotlib.pyplot as plt
from wavelets import WaveletAnalysis as WA
n = 2001
dt = 0.02
t = np.linspace(0, 10, n)
n = np.sin(-t**2) + np.random.random(n)/10
wa = WA(n, dt=dt)
scale = 1./wa.scales
power = wa.wavelet_power
T, F = np.meshgrid(t, scale)
plt.figure(figsize=(15, 6))
plt.subplot(121)
plt.plot(t, n, 'r', linewidth=1.5)
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.subplot(122)
plt.pcolormesh(T, F, power, cmap='jet')
plt.colorbar(extend='both')
plt.ylim(0.03, 1)
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Frequency [Hz]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.suptitle('Wavelet analysis', fontsize=40)
plt.show()
In [49]:
from obspy import read
import matplotlib.pyplot as plt

st = read()[: 1]
st.filter('bandpass', freqmin=0.75, freqmax=10, zerophase=True)
st.plot(color='r', grid_color='#0000AA', linewidth=2, tick_format='relative', block=False, size=(1000, 250))
Out[49]:
In [102]:
from obspy import read
import matplotlib.pyplot as plt
import numpy as np

st = read()
tr = st[2]
tr.data = tr.data - np.mean(tr.data)
tr.detrend()
tr.taper(0.005)
tr.plot()
tr.spectrogram(log=False, cmap=plt.cm.jet)
D:\Softwares\Anaconda\lib\site-packages\mkl_fft\_numpy_fft.py:158: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  output = mkl_fft.fft(a, n, axis)
In [16]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt
from obspy.signal.tf_misfit import plot_tfr

st = read()
dt = st[0].stats.delta
f1 = 1
f2 = 5

for i, tr in enumerate(st):
    st[i] = tr.filter('bandpass', freqmin=f1, freqmax=f2, zerophase=True)
    plot_tfr(st[i].data, dt=dt, fmin=f1, fmax=f2, cmap='jet')
In [24]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt
from obspy.signal.tf_misfit import plot_tf_misfits

st = read()
dt = st[0].stats.delta
f1 = 2
f2 = 8

for i, tr in enumerate(st):
    st[i] = tr.filter('bandpass', freqmin=f1, freqmax=f2, zerophase=True)

plot_tf_misfits(st[0].data, st[1].data, dt=dt, fmin=f1, fmax=f2, cmap='jet', show=False)
Out[24]:
In [13]:
from obspy.clients.syngine import Client
import numpy as np
import matplotlib.pyplot as plt

models = ['ak135f_1s',
         'ak135f_2s',
         'ak135f_5s',
         'iasp91_2s',
         'prem_i_2s',
          'prem_a_2s',
         'prem_a_5s',
          'prem_a_10s',
          'prem_a_20s'
         ]
colors = ['b', 'r', 'g', 'k', 'm', 'y', 'c']

b = 'P-50'
e = 'S+500'
client = Client()

plt.figure(figsize=(15, 6))
for i, m in enumerate(models[: 4]):
    print('model ' + str(i+1) + ': ', m )
    st = client.get_waveforms(model=m, network="IU", station="ANMO",
                          eventid="GCMT:C201002270634A", starttime=b, endtime=e)
    tr = st[0]
    #tr.filter('bandpass', freqmin=0.01, freqmax=0.15, zerophase=True)
    d = tr.data
    t = np.arange(tr.stats.npts) * tr.stats.delta
    plt.plot(t, d, linewidth=2, color=colors[i], label=m)
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.legend(fontsize=20)
plt.show()
model 1:  ak135f_1s
model 2:  ak135f_2s
model 3:  ak135f_5s
model 4:  iasp91_2s
In [9]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt

st = read()
tr = st[0]
tr.data = tr.data - np.mean(tr.data)
tr.detrend()
tr.taper(0.025)
tr.filter('bandpass', freqmin=0.5, freqmax=1.5, zerophase=True)
tr2 = tr.copy()
tr2.decimate(factor=10)
t = np.arange(len(tr.data)) * tr.stats.delta
t2 = np.arange(len(tr2.data)) * tr2.stats.delta
plt.figure(figsize=(15, 6))
plt.plot(t, tr.data, label='origin')
plt.plot(t2, tr2.data, label='downsampling')
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.legend(fontsize=20)
plt.show()
In [9]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt

st = read()
st[0].filter('bandpass', freqmin=0.5, freqmax=5, zerophase=True)
tr1 = st[0]
tr2 = st[1]
tr3 = st[2]
f11 = 1.0
f12 = 2.0
f21 = 1.2
f22 = 2.5
f31 = 2.0
f32 = 4.0
tr1.filter( 'bandpass', freqmin=f11, freqmax=f12, zerophase=True)
tr2.filter( 'bandpass', freqmin=f21, freqmax=f22, zerophase=True)
tr3.filter( 'bandpass', freqmin=f31, freqmax=f32, zerophase=True)
data1 = tr1.data
data2 = tr2.data
data3 = tr3.data
data1 = data1 / max(abs(data1))
data2 = data2 / max(abs(data2))
data3 = data3 / max(abs(data3))

t = np.linspace(0, tr1.stats.npts-1, tr1.stats.npts) * tr1.stats.delta

fig, ax = plt.subplots(figsize=(15, 9))
ax.plot( t, data1, 'g', linewidth=3, label=str(f11) + ' - ' + str(f12) + ' Hz' )
ax.text(1, 0.2, tr1.stats.channel, fontsize=20)
ax.plot( t, data2+2, 'r', linewidth=3, label=str(f21) + ' - ' + str(f22) + ' Hz' )
ax.text(1, 2.2, tr2.stats.channel, fontsize=18)
ax.plot( t, data3+4, 'b', linewidth=3, label=str(f31) + ' - ' + str(f32) + ' Hz' )
ax.text(1, 4.2, tr3.stats.channel, fontsize=18)
ax.set_xlabel( 'Time [s]', fontsize=18 )
ax.spines['top'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_yticks([])
ax.set_xticks(np.linspace(t[0], t[-1], 11))
plt.legend(fontsize=16, loc='best')
plt.xticks(fontsize=15)
plt.show()
In [72]:
from obspy import read
from obspy.signal.trigger import recursive_sta_lta, plot_trigger
import matplotlib.pyplot as plt

st = read()
tr = st.select(component="Z")[0]
tr.filter("bandpass", freqmin=1, freqmax=20)  
sta = 0.5
lta = 4
cft = recursive_sta_lta(tr.data, int(sta * tr.stats.sampling_rate),
                       int(lta * tr.stats.sampling_rate))
thrOn = 2
thrOff = 0.5
plt.figure(figsize=(15, 9))
plot_trigger(tr, cft, thrOn, thrOff)
<Figure size 1080x648 with 0 Axes>
In [56]:
from obspy import read
from obspy.signal.trigger import recursive_sta_lta, plot_trigger
import matplotlib.pyplot as plt

st = read()
tr = st.select(component="Z")[0]

tr.trigger("recstalta", sta=0.1, lta=3)
tr.plot()
Out[56]:
In [49]:
import matplotlib.pyplot as plt
from obspy.imaging.beachball import beach

np1 = [150, 87, 1]
mt = [-2.39, 1.04, 1.35, 0.57, -2.94, -0.94]

b1 = beach(np1, facecolor='b', edgecolor='k', xy=(1, 1), width=1.25)
b2 = beach(np1, facecolor='r', edgecolor='k', xy=(3, 1), width=1.5)
b3 = beach(mt, facecolor='g', edgecolor='k', xy=(7, 1), width=1.75)

plt.figure(figsize=(15, 5))
plt.plot()
ax = plt.gca()
ax.add_collection(b1)
ax.add_collection(b2)
ax.add_collection(b3)
ax.set_aspect('equal')
ax.set_xlim(0, 10)
ax.set_ylim(-1, 3)
Out[49]:
(-1, 3)
In [29]:
import matplotlib.pyplot as plt
from obspy.imaging.beachball import beach

np1 = [150, 87, 1]
mt = [-2.39, 1.04, 1.35, 0.57, -2.94, -0.94]
beach1 = beach(np1, xy=(-70, 80), width=30)
beach2 = beach(mt, xy=(50, 50), width=50)

plt.figure(figsize=(15, 8))
plt.plot([-100, 100], [0, 100], "rv", ms=80) 
ax = plt.gca()
ax.add_collection(beach1) 
ax.add_collection(beach2) 
ax.set_aspect("equal")
ax.set_xlim((-120, 120))
ax.set_ylim((-20, 120))
Out[29]:
(-20, 120)
In [17]:
from obspy.clients.iris import Client
client = Client()
result = client.traveltime(model='prem',
                           evloc=(-36.122,-72.898),
                           staloc=[(-33.45,-70.67), (47.61,-122.33), (35.69,139.69)],
                           evdepth=22.9,
                          phases = ['P', 'S', 'pP', 'sS'])
result = result.decode()
print(result)
Model: prem
Distance   Depth   Phase   Travel    Ray Param  Takeoff  Incident  Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)   (deg)    (deg)   Distance   Name 
-----------------------------------------------------------------------------------
    3.24    22.9   P         46.91    13.656     56.94    45.42     3.24   = P
    3.24    22.9   pP        51.81    13.656    123.05    45.42     3.24   = pP
    3.24    22.9   P         54.16    16.293     89.84    58.20     3.24   = P
    3.24    22.9   P         54.16    16.288     88.46    58.16     3.24   = P
    3.24    22.9   pP        57.04    16.232     94.99    57.85     3.24   = pP
    3.24    22.9   S         84.56    24.663     60.24    45.22     3.24   = S
    3.24    22.9   sS        93.17    24.664    119.75    45.22     3.24   = sS
    3.24    22.9   S         94.76    28.409     89.81    54.84     3.24   = S
    3.24    22.9   S         94.76    28.399     88.49    54.81     3.24   = S
    3.24    22.9   sS       100.42    28.308     94.83    54.55     3.24   = sS
   94.66    22.9   P        797.49     4.547     16.21    13.72    94.66   = P
   94.66    22.9   pP       804.74     4.550    163.78    13.73    94.66   = pP
   94.66    22.9   S       1469.69     8.713     17.86    14.52    94.66   = S
   94.66    22.9   sS      1482.62     8.714    162.14    14.52    94.66   = sS

In [19]:
from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 7))
ax = plot_travel_times(source_depth=10, ax=ax, fig=fig,
                       phase_list=['P', 'PP', 'S'], npoints=200)
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
In [62]:
from obspy.taup import TauPyModel
import matplotlib.pyplot as plt

phase_list=['PP', 'SSS', 'SKS', 'PKP', 'SKKS', 'PKKP', 'SKIKS', 'PKJKP', 'SKJKS', 'PKIKP', 'PKIIKP']
model = TauPyModel(model='iasp91')
arrivals = model.get_ray_paths(500, 140, phase_list=phase_list)
fig, ax = plt.subplots(figsize=(12, 8))
ax.xaxis.set_ticks_position('top')
ax.invert_yaxis()

ax = arrivals.plot_rays(plot_type='cartesian', phase_list=phase_list,
                   plot_all=False, legend=True, ax = ax)
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
In [38]:
from obspy.taup import TauPyModel
import matplotlib.pyplot as plt

model = TauPyModel(model='iasp91')
arrivals = model.get_ray_paths(500, 140, phase_list=['Pdiff', 'SS'])
fig, ax = plt.subplots(figsize=(12, 10))
arrivals.plot_rays(plot_type='spherical', phase_list=['Pdiff', 'SS'],
                   legend=True, fig=fig)
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
Out[38]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0x1fbeb58de10>
In [39]:
from obspy.taup.tau import plot_ray_paths
import matplotlib.pyplot as plt

fig, ax = plt.subplots(subplot_kw=dict(polar=True), figsize=(10, 9))
ax = plot_ray_paths(source_depth=100, ax=ax, fig=fig, phase_list=['P', 'PKP'],
                    npoints=25)
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
In [41]:
import numpy as np
import matplotlib.pyplot as plt

from obspy.taup import TauPyModel


PHASES = [
    # Phase, distance
    ('P', 26),
    ('PP', 60),
    ('PPP', 94),
    ('PPS', 155),
    ('p', 3),
    ('pPcP', 100),
    ('PKIKP', 170),
    ('PKJKP', 194),
    ('S', 65),
    ('SP', 85),
    ('SS', 134.5),
    ('SSS', 204),
    ('p', -10),
    ('pP', -37.5),
    ('s', -3),
    ('sP', -49),
    ('ScS', -44),
    ('SKS', -82),
    ('SKKS', -120),
]

model = TauPyModel(model='iasp91')

fig, ax = plt.subplots(subplot_kw=dict(polar=True), figsize=(10, 10))

# Plot all pre-determined phases
for phase, distance in PHASES:
    arrivals = model.get_ray_paths(700, distance, phase_list=[phase])
    ax = arrivals.plot_rays(plot_type='spherical',
                            legend=False, label_arrivals=True,
                            plot_all=True,
                            show=False, ax=ax)

# Annotate regions
ax.text(0, 0, 'Solid\ninner\ncore',
        horizontalalignment='center', verticalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
ocr = (model.model.radius_of_planet -
       (model.model.s_mod.v_mod.iocb_depth +
        model.model.s_mod.v_mod.cmb_depth) / 2)
ax.text(np.deg2rad(180), ocr, 'Fluid outer core',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
mr = model.model.radius_of_planet - model.model.s_mod.v_mod.cmb_depth / 2
ax.text(np.deg2rad(180), mr, 'Solid mantle',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))

plt.show()
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
In [42]:
from __future__ import print_function

import obspy
import obspy.clients.fdsn


client = obspy.clients.fdsn.Client("EMSC")

events = client.get_events(minlatitude=46.1, maxlatitude=46.3,
                           minlongitude=7.6, maxlongitude=7.8,
                           starttime=obspy.UTCDateTime("2012-04-03"),
                           endtime=obspy.UTCDateTime("2012-04-04"))

print("found %s event(s):" % len(events))
for event in events:
    print(event)
found 2 event(s):
Event:	2012-04-03T02:47:32.500000Z | +46.220,   +7.660 | 1.8 ml | manual

	            resource_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000005")
	          creation_info: CreationInfo(agency_uri=ResourceIdentifier(id="smi:smi-registry/organization/EMSC"), author_uri=ResourceIdentifier(id="smi:smi-registry/organization/ZUR"), creation_time=UTCDateTime(2014, 12, 1, 12, 32, 21, 193214))
	    preferred_origin_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000005/origin/LAB57")
	 preferred_magnitude_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000005/origin/LAB57/mag/1")
	                   ---------
	     event_descriptions: 1 Elements
	                origins: 1 Elements
	             magnitudes: 1 Elements
Event:	2012-04-03T02:45:03.300000Z | +46.220,   +7.710 | 2.5 ml | automatic

	            resource_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000004")
	          creation_info: CreationInfo(agency_uri=ResourceIdentifier(id="smi:smi-registry/organization/EMSC"), author_uri=ResourceIdentifier(id="smi:smi-registry/organization/ZUR"), creation_time=UTCDateTime(2014, 12, 1, 14, 41, 41, 514185))
	    preferred_origin_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000004/origin/ee086")
	 preferred_magnitude_id: ResourceIdentifier(id="quakeml:eu.emsc/event/20120403_0000004/origin/ee086/mag/1")
	                   ---------
	     event_descriptions: 1 Elements
	                origins: 1 Elements
	             magnitudes: 1 Elements
In [14]:
from obspy.taup.taup_time import TauPTime as TPT

model = TPT(model='prem', phase_list=['P', 'S'], depth=50, degrees=60, receiver_depth=0)

print(model.source_depth)
50
In [ ]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt

st = read('1*.sac')
tr = st[3]
dt = tr.stats.delta
fs = 1 / dt
n = tr.stats.npts

print(dt)
tr.detrend()
tr.data = tr.data - np.mean(tr.data)
tr.taper(0.025)
tr.filter('bandpass', freqmin=0.05, freqmax=0.1, corners=4, zerophase=True)
t = np.arange(n) * dt + tr.stats.sac.b
s = tr.data
s = s /max(abs(s))


fft_s = abs(np.fft.fft(s))
f = np.arange(n) * fs / n

fft_s = fft_s[0: n//2]
f = f[0: n//2]

plt.figure(figsize=(20, 10))
plt.subplot(211)
plt.plot(t, s, 'k', linewidth=2)
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

plt.subplot(212)
plt.stem(f, fft_s, 'b', linewidth=2)
plt.xlabel('Frequency [Hz]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()
0.04
D:\Softwares\Anaconda\lib\site-packages\matplotlib\cbook\deprecation.py:107: MatplotlibDeprecationWarning: stem() got an unexpected keyword argument 'linewidth'. This will raise a TypeError in future versions.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [5]:
from obspy import read
import numpy as np
import matplotlib.pyplot as plt
from obspy.signal.filter import envelope

st = read('STA*.SAC')
tr = st[0]
n0 = 100000
tr.data = tr.data[: n0]
tr.detrend()
tr.data = tr.data - np.mean(tr.data)
tr.taper(0.05)
f1 = 0.2
f2 = 1.5
tr.filter('bandpass', freqmin=f1, freqmax=f2, zerophase=True)
s = tr.data
tr.filter('bandpass', freqmin=f1, freqmax=f2*0.75, zerophase=True)
s2 = tr.data

s3 = envelope(s2)

t = np.arange(len(s)) * tr.stats.delta + tr.stats.sac.b

f = np.arange(n0) / tr.stats.delta / n0
fft_s = abs(np.fft.fft(s))
f = f[0: n0//2]
fft_s = fft_s[0: n0//2]

plt.figure(figsize=(18, 10))
plt.subplot(211)
plt.plot(t, s, 'r', linewidth=1.5, label='filtered 1')
plt.plot(t, s2, 'b', linewidth=1.5, label='filtered 2')
plt.plot(t, s3, 'g', linewidth=3, label='envelope')
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.legend(fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

plt.subplot(212)
plt.plot(f, fft_s, 'b', linewidth=1.5)
plt.xlabel('Frequency [Hz]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xlim(f1-f1*0.5, f2+f1*0.5)
plt.show()
D:\Softwares\Anaconda\lib\site-packages\obspy\signal\headers.py:93: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  ], align=True)
D:\Softwares\Anaconda\lib\site-packages\IPython\core\pylabtools.py:125: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
In [30]:
import numpy as np
from obspy import read
from obspy.signal import PPSD
import matplotlib.pyplot as plt


st = read('1*.sac')
tr = st[0]
paz = {'gain': 60077000.0,
       'poles': [-0.037004+0.037016j, -0.037004-0.037016j,
                 -251.33+0j, -131.04-467.29j, -131.04+467.29j],
       'sensitivity': 2516778400.0,
       'zeros': [0j, 0j]}
ppsd = PPSD(tr.stats, paz)
ppsd.add(st) 

plt.figure(figsize=(15, 10))
ppsd.plot(cmap=plt.cm.jet) 
<Figure size 1080x720 with 0 Axes>
In [35]:
from obspy import  read
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d as interp2d

st = read('Australis_data\Data\*BHZ*.SAC')
nr = 20
st = st[: nr]
t1 = 800
t2 = 3000
dt = st[0].stats.delta
n1 = int(t1/dt)
n2 = int(t2/dt)
tt = np.arange(n2-n1) * dt + t1

plt.figure(figsize=(20, 10))
for i, tr in enumerate(st):
    dt = tr.stats.delta
    d = tr.data[n1: n2]
    tr.data = d - np.mean(d)
    tr.detrend()
    tr.taper(0.0025)
    tr.filter('bandpass', freqmin=0.005, freqmax=0.05, corners=4, zerophase=True)
    tr.data = tr.data / max(abs(tr.data)) * 1.75
    n = len(tr.data)
    t = np.arange(n) * dt + t1
    d = tr.data
    d = d / max(abs(d)) * 5
    plt.plot(t, d+tr.stats.sac.gcarc, color='r', alpha=0.25, linewidth=2)

plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Distance [$^o$]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Travel times of some phases
phase1 = dict()
phase2 = dict()
colors = ['#FF0000', '#00FF00', '#0000FF', '#00FFFF', '#000000', '#FFFF00', '#FF00FF', '#999999', '#FFA500', '#800080']
for i, tr in enumerate(st):
    n = len(tr.data)
    t = np.arange(n) * dt + t1
    phase1['0'] = tr.stats.sac.t0
    phase1['1'] = tr.stats.sac.t1
    phase1['2'] = tr.stats.sac.t2
    phase1['3'] = tr.stats.sac.t3
    phase1['4'] = tr.stats.sac.t4
    phase1['5'] = tr.stats.sac.t5
    phase1['6'] = tr.stats.sac.t6
    phase1['7'] = tr.stats.sac.t7
    phase1['8'] = tr.stats.sac.t8
    phase1['9'] = tr.stats.sac.t9
    phase2['0'] = tr.stats.sac.kt0
    phase2['1'] = tr.stats.sac.kt1
    phase2['2'] = tr.stats.sac.kt2
    phase2['3'] = tr.stats.sac.kt3
    phase2['4'] = tr.stats.sac.kt4
    phase2['5'] = tr.stats.sac.kt5
    phase2['6'] = tr.stats.sac.kt6
    phase2['7'] = tr.stats.sac.kt7
    phase2['8'] = tr.stats.sac.kt8
    phase2['9'] = tr.stats.sac.kt9
    if i < 1:
        index = 0
        for (k1, k2) in zip(phase1.keys(), phase2.keys()):
            plt.plot(phase1[k1], tr.stats.sac.gcarc, 'o', color=colors[index], linewidth=1, label=phase2[k2])
            index += 1
    else:
        index = 0
        for k in phase1.keys():
            plt.plot(phase1[k], tr.stats.sac.gcarc, 'o', color=colors[index],  linewidth=1)
            index += 1
plt.legend(loc='best', fontsize=13)
plt.xlim(t1, t2)

plt.show()
In [315]:
from obspy import  read
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d as interp2d

st = read('Australis_data\Data\*BHZ*.SAC')
nr = 20
st = st[: nr]
t1 = 800
t2 = 1600
dt = st[0].stats.delta
n1 = int(t1/dt)
n2 = int(t2/dt)
tt = np.arange(n2-n1) * dt + t1
seis = np.zeros((nr, (n2-n1)))
dist = []

for i, tr in enumerate(st):
    dt = tr.stats.delta
    dist.append(tr.stats.sac.gcarc)
    d = tr.data[n1: n2]
    tr.data = d - np.mean(d)
    tr.detrend()
    tr.taper(0.0025)
    tr.filter('bandpass', freqmin=0.005, freqmax=0.05, corners=4, zerophase=True)
    tr.data = tr.data / max(abs(tr.data)) * 1.75
    n = len(tr.data)
    t = np.arange(n) * dt + t1
    seis[i] = tr.data

plt.figure(figsize=(15, 12))
T, D = np.meshgrid(tt, dist)
plt.scatter(T, D, c=seis, marker='s', s=500, cmap='RdGy_r', edgecolor='None')
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Distance [$^o$]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
for i, tr in enumerate(st):
    n = len(tr.data)
    t = np.arange(n) * dt + t1
    #plt.plot(t, seis[i]+tr.stats.sac.gcarc, color='#ECCBAF', linewidth=2)
    plt.plot(t, seis[i]+tr.stats.sac.gcarc, color='g', linewidth=2)
    if i < 1:
        plt.plot(tr.stats.sac.t3, tr.stats.sac.gcarc, 'ro', linewidth=1, label=tr.stats.sac.kt3)
        plt.plot(tr.stats.sac.t4, tr.stats.sac.gcarc, 'bo', linewidth=1, label=tr.stats.sac.kt4)
        plt.plot(tr.stats.sac.t5, tr.stats.sac.gcarc, 'go', linewidth=1, label=tr.stats.sac.kt5)
    else:
        plt.plot(tr.stats.sac.t3, tr.stats.sac.gcarc, 'ro', linewidth=1)
        plt.plot(tr.stats.sac.t4, tr.stats.sac.gcarc, 'bo', linewidth=1)
        plt.plot(tr.stats.sac.t5, tr.stats.sac.gcarc, 'go', linewidth=1)
plt.legend(loc='best', fontsize=13)
plt.xlim(t1, t2)

plt.show()
In [279]:
from obspy.taup import TauPyModel as TP
import numpy as np
import matplotlib.pyplot as plt

n1 = 25
n2 = 25
times = np.zeros((n1, n2))
dist = np.arange(n1) * 1 + 130
depth = np.arange(n2) * 2 + 20
tp = TP('prem')
for i, d1 in enumerate(dist):
    for j, d2 in enumerate(depth):
        arr = tp.get_travel_times(source_depth_in_km=d2, distance_in_degree=d1, phase_list=['PKIKP', 'PKIIKP'])
        times[i, j] = arr[1].time
        #print(d1, d2, times[i, j])

plt.figure(figsize=(15, 12))
plt.subplot(211)
#plt.plot(times, dist, 'ro', linewidth=2)
#plt.xlabel('Travel times [s]', fontsize=25)
#plt.ylabel('$ Distance [^o]$', fontsize=25)
plt.pcolormesh(depth, dist, times, cmap='RdGy')

plt.ylabel('Distance [deg]', fontsize=25)
plt.xlabel('Depth [km]', fontsize=25)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.colorbar()

plt.subplot(212)
plt.plot(depth, times[0], 'ko', linewidth=2)
plt.ylabel('Time [s]', fontsize=25)
plt.xlabel('Depth [km]', fontsize=25)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()
D:\Softwares\Anaconda\lib\site-packages\obspy\taup\tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
In [11]:
import numpy as np
import matplotlib.pyplot as plt

from obspy.imaging.cm import obspy_sequential
from obspy.signal.array_analysis import array_transff_wavenumber

# generate array coordinates
r = 300
coords1 = np.array([[0., r, 0.], 
                   [r*np.sqrt(3)/2, r/2., 0.],
                   [r*np.sqrt(3)/2, -r/2., 0.],
                   [0., -r, 0.],
                   [-r*np.sqrt(3)/2, -r/2, 0.],
                   [-r*np.sqrt(3)/2, r/2, 0],
                  ])
coords2 = r * np.array([[-3, 0, 0],
                    [-1.5, 0, 0],
                    [0, 0, 0],
                    [1.5, 0, 0],
                    [3, 0, 0],
                    [0, -3, 0],
                    [0, -1.5, 0],
                    [0, 1.5, 0],
                    [0, 3, 0]]) / 3
n3 = 10
theta = np.linspace(-180, 180, n3) / 180 * np.pi
coords3 = np.zeros((n3, 3))
coords3[:, 0] = np.cos(theta) * r
coords3[:, 1] = np.sin(theta) * r
n4 = 21
coords4 = np.zeros((n4, 3))
coords4[:, 0] = np.linspace(-r, r, n4)
coords4[:, 1] = np.cos(coords4[:, 0]/r*4) * r
# coordinates in km
coords1 /= 1000
coords2 /= 1000
coords3 /= 1000
coords4 /= 1000

# set limits for wavenumber differences to analyze
klim = 100.
kxmin = -klim
kxmax = klim
kymin = -klim
kymax = klim
kstep = klim / 50.

# compute transfer function as a function of wavenumber difference
transff1 = array_transff_wavenumber(coords1, klim, kstep, coordsys='xy')
transff2 = array_transff_wavenumber(coords2, klim, kstep, coordsys='xy')
transff3 = array_transff_wavenumber(coords3, klim, kstep, coordsys='xy')
transff4 = array_transff_wavenumber(coords4, klim, kstep, coordsys='xy')

# Plot array1
plt.figure(figsize=(24, 10))
# Array map
plt.subplot(241)
plt.scatter(coords1[:, 0], coords1[:, 1], marker='v', s=200, facecolor='b')
plt.xlabel('X [km]', fontsize=15)
plt.ylabel('Y [km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.subplot(242)
plt.scatter(coords2[:, 0], coords2[:, 1], marker='v', s=200, facecolor='b')
plt.xlabel('X [km]', fontsize=15)
plt.ylabel('Y [km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.subplot(243)
plt.scatter(coords3[:, 0], coords3[:, 1], marker='v', s=200, facecolor='b')
plt.xlabel('X [km]', fontsize=15)
plt.ylabel('Y [km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.subplot(244)
plt.scatter(coords4[:, 0], coords4[:, 1], marker='v', s=200, facecolor='b')
plt.xlabel('X [km]', fontsize=15)
plt.ylabel('Y [km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Array Transfer Function
plt.subplot(245)
plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
           transff1.T, cmap='jet')
plt.colorbar()
plt.xlabel('Slowness [s/km]', fontsize=15)
plt.ylabel('Slowness [s/km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Plot array2
plt.subplot(246)
plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
           transff2.T, cmap='viridis')
plt.colorbar()
plt.xlabel('Slowness [s/km]', fontsize=15)
plt.ylabel('Slowness [s/km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Plot array3
plt.subplot(247)
plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
           transff3.T, cmap='hot')
plt.colorbar()
plt.xlabel('Slowness [s/km]', fontsize=15)
plt.ylabel('Slowness [s/km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Plot array4
plt.subplot(248)
plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,
           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,
           transff4.T, cmap='bone')
plt.colorbar()
plt.xlabel('Slowness [s/km]', fontsize=15)
plt.ylabel('Slowness [s/km]', fontsize=15)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.suptitle('Array Transfer Function', fontsize=25)
plt.show()
In [ ]:
from obspy.clients.syngine import Client

client = Client()

st = client.get_waveforms(model="ak135f_5s", network="IU", station="ANMO",
                          eventid="GCMT:C201002270634A")
print(st)  

st[0].plot() 
In [90]:
from obspy.clients.fdsn import Client
from obspy import UTCDateTime

b = UTCDateTime('2019-08-05T00:00:00.000')
e = UTCDateTime('2599-12-31T23:59:59.999')

clint = Client('IRIS')
sta = clint.get_stations(minlongitude=120,
                         maxlongitude=130,
                         minlatitude=20,
                         maxlatitude=50,
                       station='*',
                       starttime=b,
                       endtime=e)
print(sta)
stations = sta.get_contents()['stations']
for s in stations[: 10]:
    sn = s.split()[0]
    print(sn)
Inventory created at 2019-08-11T09:23:15.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?starttime=2019-08-05...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (9):
			CB, IC, IM, IU, JP, KG, KS, SY, TW
		Stations (106):
			CB.CN2 (Changchun,Jilin Province)
			CB.DL2 (Dalian,Liaoning Province)
			CB.SNY (Shenyang,Liaoning Province)
			IC.MDJ (Mudanjiang, Heilongjiang Province, China)
			IC.SSE (Shanghai, China)
			IM.KSAR (Wonju Array Beam, South Korea)
			IU.INCN (Inchon, Republic of Korea)
			IU.TATO (Taipei, Taiwan)
			JP.JMJ (Miyakojima Island)
			JP.JMJ2 (Miyakojima Island)
			JP.JOW (Okinawa Kunigami)
			JP.JTU (Tsushima Kamiagata)
			JP.YOJ (Yonagunijima Island)
			KG.TJN (Taejon, Republic of Korea)
			KS.BUS2 (Busan, South Korea)
			KS.CHJ2 (Chungju, South Korea)
			KS.NAWB (Namwon, South Korea)
			KS.SEHB (Seawha, South Korea)
			KS.SEO2 (Seoul, South Korea)
			SY.ANPB (ANPB synthetic)
			SY.BUS2 (BUS2 synthetic)
			SY.CHGB (CHGB synthetic)
			SY.CHJ2 (CHJ2 synthetic)
			SY.CN2 (CN2 synthetic)
			SY.DL2 (DL2 synthetic)
			SY.HGSD (HGSD synthetic)
			SY.HWAB (HWAB synthetic)
			SY.INCN (INCN synthetic)
			SY.JGPD (JGPD synthetic)
			SY.JOW (JOW synthetic)
			SY.JTU (JTU synthetic)
			SY.LYUB (LYUB synthetic)
			SY.MDJ (MDJ synthetic)
			SY.MDPD (MDPD synthetic)
			SY.NACB (NACB synthetic)
			SY.NAWB (NAWB synthetic)
			SY.NE00 (NE00 synthetic)
			SY.NE01 (NE01 synthetic)
			SY.NE02 (NE02 synthetic)
			SY.NE03 (NE03 synthetic)
			SY.NE04 (NE04 synthetic)
			SY.NE05 (NE05 synthetic)
			SY.NE06 (NE06 synthetic)
			SY.NE07 (NE07 synthetic)
			SY.NE08 (NE08 synthetic)
			SY.NE09 (NE09 synthetic)
			SY.NE10 (NE10 synthetic)
			SY.NE11 (NE11 synthetic)
			SY.NE12 (NE12 synthetic)
			SY.NE14 (NE14 synthetic)
			SY.NE15 (NE15 synthetic)
			SY.NE16 (NE16 synthetic)
			SY.NE17 (NE17 synthetic)
			SY.NE18 (NE18 synthetic)
			SY.NE19 (NE19 synthetic)
			SY.NE20 (NE20 synthetic)
			SY.NE21 (NE21 synthetic)
			SY.NE22 (NE22 synthetic)
			SY.NE23 (NE23 synthetic)
			SY.NE24 (NE24 synthetic)
			SY.NE25 (NE25 synthetic)
			SY.NE26 (NE26 synthetic)
			SY.NE27 (NE27 synthetic)
			SY.NE28 (NE28 synthetic)
			SY.NE29 (NE29 synthetic)
			SY.NE30 (NE30 synthetic)
			SY.NE31 (NE31 synthetic)
			SY.NE32 (NE32 synthetic)
			SY.NE33 (NE33 synthetic)
			SY.NE34 (NE34 synthetic)
			SY.NE35 (NE35 synthetic)
			SY.NE36 (NE36 synthetic)
			SY.PDBD (PDBD synthetic)
			SY.PSRD (PSRD synthetic)
			SY.RLNB (RLNB synthetic)
			SY.SEHB (SEHB synthetic)
			SY.SEO2 (SEO2 synthetic)
			SY.SHRD (SHRD synthetic)
			SY.SMSD (SMSD synthetic)
			SY.SNY (SNY synthetic)
			SY.SSE (SSE synthetic)
			SY.SSLB (SSLB synthetic)
			SY.TATO (TATO synthetic)
			SY.TDCB (TDCB synthetic)
			SY.TJN (TJN synthetic)
			SY.TPUB (TPUB synthetic)
			SY.TWGB (TWGB synthetic)
			SY.TWKB (TWKB synthetic)
			SY.WFSB (WFSB synthetic)
			SY.YHNB (YHNB synthetic)
			SY.YOJ (YOJ synthetic)
			SY.YULB (YULB synthetic)
			TW.CHGB (Ching-ging)
			TW.HGSD (Houtzshan)
			TW.LYUB (Lanyu)
			TW.NACB (Ninganchiao)
			TW.RLNB (Erlin)
			TW.SSLB (Suanglung)
			TW.TATO (TAI-PEI)
			TW.TDCB (Techi)
			TW.TPUB (Dapu)
			TW.TWGB (Taitung)
			TW.TWKB (Kenting)
			TW.WFSB (Wufenshan)
			TW.YHNB (Yeheng)
			TW.YULB (Yuli)
		Channels (0):

CB.CN2
CB.DL2
CB.SNY
IC.MDJ
IC.SSE
IM.KSAR
IU.INCN
IU.TATO
JP.JMJ
JP.JMJ2
In [309]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
nf = 20
shift = np.random.uniform(-1, 1, nf) * np.pi
ss = np.zeros((nf, 1001))
for i in range(nf):
    ss[i] = np.sin((i+1)*2*np.pi*t+shift[i]) * np.random.rand(1001)
s = sum(ss)
#s = s1

fft_s = abs(np.fft.fft(s))
f = np.arange(1001) * 100 / 1001
fft_s = fft_s[: 500]
f = f[: 500]

plt.figure(figsize=(12, 10))

plt.subplot(311)
plt.plot(t, s, linewidth=2)
plt.subplot(312)
plt.plot(f, fft_s)
plt.subplot(313)
plt.plot(np.arange(nf), shift)
plt.plot(np.arange(nf), shift, 'ko')
plt.show()
In [19]:
import numpy as np
import matplotlib.pyplot as plt

n = 21
x = np.linspace(0, 50, n)
y = np.linspace(0, 50, n)
X, Y = np.meshgrid(x, y)
s = np.sin(X*X)*np.cos(Y) + np.cos(Y)*np.sin(X) + np.random.random((n, n)) + np.random.random((n, n)) + \
    np.random.random((n, n)) +np.random.random((n, n)) + np.random.random((n, n)) + np.random.random((n, n))
plt.figure(figsize=(15,10))
plt.scatter(X, Y, c=s, marker='s',linewidth=1, s=300, alpha=1, cmap='viridis', edgecolor='None')
plt.xlim(-5,55)
plt.ylim(-5,55)
#plt.contourf(X, Y, s, cmap='viridis')
plt.colorbar()
plt.show()
In [22]:
import numpy as np
import matplotlib.pyplot as plt

n = 21
x = np.linspace(0, 50, n)
y = np.linspace(0, 50, n)
X, Y = np.meshgrid(x, y)
s = np.sin(X*X)*np.cos(Y) + np.cos(Y)*np.sin(X) + np.random.random((n, n)) + np.random.random((n, n)) +\
np.random.random((n, n)) +np.random.random((n, n)) + np.random.random((n, n)) + np.random.random((n, n))
#plt.scatter(X, Y, c=s, marker='s',linewidth=1, s=120, alpha=1, cmap='viridis', edgecolor='None')
plt.figure(figsize=(15,10))
plt.xlim(-5,55)
plt.ylim(-5,55)
plt.contourf(X, Y, s, cmap='viridis')
plt.colorbar()
plt.title('contourf', fontsize=20)
plt.show()
In [24]:
import numpy as np
import matplotlib.pyplot as plt

arr = np.random.random((10, 10))
plt.figure(figsize=(18,10))
fig, axs = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
for ax in axs.flatten():
    im = ax.pcolormesh(arr,cmap='viridis')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
fig.colorbar(im, ax=axs, shrink=0.75, extend='both')
plt.suptitle('plt.pcolormesh', fontsize=25)
plt.show()
<Figure size 1296x720 with 0 Axes>
In [25]:
import numpy as np
import matplotlib.pyplot as plt
arr = np.random.random((100, 100))
x = np.linspace(-5,5, 51)
y = np.linspace(-3, 3, 61)
X, Y = np.meshgrid(x, y)
v = np.exp(-X**2-Y**2)
r = np.sqrt(X**2+Y**2)
plt.figure(figsize=(20, 5))
plt.subplot(131)
plt.pcolormesh(arr, cmap='RdBu_r')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.colorbar()
plt.subplot(132)
plt.pcolormesh(X, Y, v, cmap='RdBu_r')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.colorbar(extend='both')
plt.subplot(133)
plt.pcolormesh(X, Y, r, cmap='seismic')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.colorbar(extend='both')
plt.suptitle('pcolormesh different colorbars', fontsize=30)
plt.show()
In [32]:
import numpy as np
import matplotlib.pyplot as plt

n = 41
x = np.linspace(-10, 10, n)
y = np.linspace(-10, 10, n)
X, Y = np.meshgrid(x, y)
v = np.sqrt(X**2+Y**2)
plt.figure(figsize=(10,8))
plt.pcolormesh(X, Y, v, cmap='hot')
plt.xlabel('X axis', fontsize=20)
plt.ylabel('Y axis', fontsize=20)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
cbar = plt.colorbar(extend='both', shrink=0.75)
cbar.set_label('Power', fontsize=20)
cbar.set_ticks(np.arange(0, 15, 2))
plt.title('Set label of colorbar', fontsize=20)
plt.show()
In [36]:
import numpy as np
import matplotlib.pyplot as plt

n = 21
x = np.linspace(0, 50, n)
y = np.linspace(0, 50, n)
X, Y = np.meshgrid(x, y)
s = np.sin(X*X)*np.cos(Y) + np.cos(Y)*np.sin(X) + np.random.random((n, n)) + np.random.random((n, n)) + \
np.random.random((n, n)) +np.random.random((n, n)) + np.random.random((n, n)) + np.random.random((n, n))
plt.figure(figsize=(10, 8))
plt.scatter(X, Y, c=s, marker='o',linewidth=1, s=abs(s)*100, alpha=1, cmap='seismic', edgecolor='None')
plt.xlim(-5,55)
plt.ylim(-5,55)
#plt.contourf(X, Y, s, cmap='viridis')
plt.colorbar(extend='both', shrink=0.75)
plt.title('different sizes of scatters', fontsize=20)
plt.show()
In [38]:
import numpy as np
import matplotlib.pyplot as plt
n = 6
x1 = np.linspace(0, 50, n)
x2 = np.linspace(1, 51, n)
y1 = np.linspace(0, 50, n)
y2 = np.linspace(1, 51, n)
X1, Y1 = np.meshgrid(x1, y1)
X2, Y2 = np.meshgrid(x2, y2)
s1 = np.sin(X1*X1)*np.cos(Y1) + np.cos(Y1)*np.sin(X1) + np.random.random((n,n))
s2 = np.random.random((n, n))
plt.figure(figsize=(10, 8))
plt.subplot(111)
plt.scatter(X1, Y1, c=s1, marker='s',linewidth=1, s=500, alpha=1, cmap='viridis', edgecolor='None', label='squares')
plt.colorbar(extend='both', shrink=0.75)
plt.subplot(111)
plt.scatter(X2, Y2, c=s2, marker='o',linewidth=1, s=500, alpha=1, cmap='hot', edgecolor='None',label='circles')
plt.xlim(-5,55)
plt.ylim(-5,55)
plt.colorbar(extend='both', shrink=0.75)
plt.legend(loc='best')
plt.show()
D:\Anaconda\lib\site-packages\matplotlib\cbook\deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [41]:
import numpy as np
import matplotlib.pyplot as plt

na = 37
nr = 7
azimuths = np.radians(np.linspace(0, 360, na))
zeniths = np.linspace(0, 10, nr)
r, theta = np.meshgrid(zeniths, azimuths)
values = np.random.random((na, nr))
fig, ax = plt.subplots(subplot_kw=dict(projection='polar'), figsize=(10, 8))
p = ax.contourf(theta, r, values, cmap='plasma')
ax.set_title('Polar contourf plot', fontsize=25)
cbar = plt.colorbar(p, extend='both', shrink=0.75)
plt.show()
In [45]:
import numpy as np
import matplotlib.pyplot as plt

azimuths = np.radians(np.linspace(0, 360, 37))
#zeniths = np.arange(0, 11, 10)
zeniths = np.linspace(0, 10, 2)
r, theta = np.meshgrid(zeniths, azimuths)
values1 = np.random.random((azimuths.size, zeniths.size))
values2 = np.random.random((azimuths.size, zeniths.size))
fig, axs = plt.subplots(1, 2, subplot_kw=dict(projection='polar'), figsize=(12, 9))
p1 = axs[0].contourf(theta, r, values1, 100, cmap='jet')
plt.colorbar(p1, ax=axs[1])
p2 = axs[1].contourf(theta, r, values2, 100, cmap='viridis')
plt.colorbar(p2, ax=axs[1])
plt.suptitle('Polar contourf', fontsize=20)
plt.show()
In [87]:
import numpy as np
import matplotlib.pyplot as plt

theta = np.arange(0, 2*np.pi, np.pi/360)
plt.figure(figsize=(15, 6))
plt.subplot(121, polar=True)
plt.plot(theta, 1.4*np.cos(5*theta), 'b--', linewidth=2)
plt.plot(theta, 1.8*np.cos(4*theta), 'g', linewidth=2)
plt.rgrids(np.linspace(0.5, 2.0, 4), angle=45, color='r', fontsize=15)
plt.thetagrids([0, 90, 180, 270], fontsize=15)
plt.subplot(122, polar=True)
plt.plot(theta, 1.6*np.ones_like(theta), 'b', linewidth=2) #绘制同心圆
plt.plot(3*theta, theta/3, 'g--', linewidth=2)
plt.rgrids(np.linspace(0.5, 2.0, 4), angle=45, color='r', fontsize=15)
plt.thetagrids([0, 90, 180, 270], fontsize=15)
plt.show()
In [337]:
import numpy as np
import matplotlib.pyplot as plt

n = 360
theta = np.arange(n) * np.pi / 18
rho = np.random.rand(n) + 0.5
area = 300 * rho ** 2

plt.figure(figsize=(20, 7))
plt.subplot(131, projection='polar')
plt.scatter(theta, rho, c=theta, s=area, cmap='jet', alpha=0.75)
plt.title('poloar 1')

ax = plt.subplot(132, projection='polar')
ax.scatter(theta, rho, c=theta, s=area, cmap='RdGy_r', alpha=0.75)
ax.set_thetamin(90)
ax.set_thetamax(180)
plt.title('poloar 2')

ax = plt.subplot(133, projection='polar')
ax.scatter(theta, rho, c=theta, s=area, cmap='viridis', alpha=0.75)
ax.set_rorigin(-0.75)
ax.set_theta_zero_location('W', offset=10)
plt.title('poloar 3')
plt.show()
In [349]:
import numpy as np
import matplotlib.pyplot as plt

n = 360
theta = np.arange(n) * np.pi / 18
rho = np.random.rand(n) + 0.5
area = 300 * rho ** 2

plt.figure(figsize=(20, 7))
plt.subplot(141, projection='polar')
plt.scatter(theta, rho, c=theta, s=area, cmap='hot', alpha=0.75, marker=r'$\clubsuit$')
plt.title('poloar 1')

ax = plt.subplot(142, projection='polar')
ax.scatter(theta, rho, c=theta, s=area, cmap='jet', alpha=0.75, marker=r'$\spadesuit$')
ax.set_thetamin(90)
ax.set_thetamax(180)
plt.title('poloar 2')

ax = plt.subplot(143, projection='polar')
ax.scatter(theta, rho, c=theta, s=area, cmap='plasma', alpha=0.75, marker=r'$\heartsuit$')
ax.set_rorigin(-0.75)
ax.set_theta_zero_location('W', offset=10)
plt.title('poloar 3')

ax = plt.subplot(144, projection='polar')
ax.scatter(theta, rho, c=theta, s=area, cmap='Spectral', alpha=0.75, marker=r'$\diamondsuit$')
ax.set_rorigin(-0.75)
ax.set_theta_zero_location('W', offset=10)
plt.title('poloar 3')
plt.show()
In [83]:
import numpy as np
import matplotlib.pyplot as plt

n = 8
X, Y = np.mgrid[0:n, 0:n]
T = np.arctan2(Y - n / 2., X - n/2.)
R = 10 + np.sqrt((Y - n / 2.0) ** 2 + (X - n / 2.0) ** 2)
U, V = R * np.cos(T), R * np.sin(T)
plt.figure(figsize=(22, 6))
plt.subplot(131)
plt.quiver(X, Y, U, V, R, alpha=.5, cmap='jet')
plt.xlim(-1, n)
plt.xticks(())
plt.ylim(-1, n)
plt.yticks(())

plt.subplot(132)
plt.quiver(X, Y, U, V, R, edgecolor='k', linewidth=2, cmap='viridis')
plt.xlim(-1, n)
plt.xticks(())
plt.ylim(-1, n)
plt.yticks(())

plt.subplot(133)
plt.quiver(X, Y, U*10, V, edgecolor='k', facecolor='None', linewidth=1)
plt.xlim(-1, n)
plt.xticks(())
plt.ylim(-1, n)
plt.yticks(())
plt.suptitle('Flow arrows', fontsize=28)
plt.show()
In [89]:
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(10, 6))
ax = plt.subplot(111)
t = np.linspace(0, 5, 501)
s = np.cos(2*np.pi*t)
line, = plt.plot(t, s,'k', lw=4)
plt.annotate('local max', xy=(2, 1), xytext=(4, 1.5), arrowprops=dict(edgecolor='b', facecolor=
'red', shrink=0.05), fontsize=15)
plt.ylim(-2, 2)
plt.title('Arrow marks', fontsize=25)
plt.show()
In [95]:
import numpy as np
import matplotlib.pyplot as plt

with plt.style.context(('dark_background')):
    plt.figure(figsize=(10,6))
    plt.plot(np.sin(np.linspace(0, 4*np.pi, 101)**2)*np.exp(-(np.linspace(0, 4*np.pi, 101))), 'r-^')
    plt.title('black background', fontsize=25)
    plt.grid(True)
plt.show()
In [354]:
import numpy as np
import matplotlib.pyplot as plt


category_names = ['Strongly disagree', 'Disagree',
                  'Neither agree nor disagree', 'Agree', 'Strongly agree']
results = {
    'Question 1': [10, 15, 17, 32, 26],
    'Question 2': [26, 22, 29, 10, 13],
    'Question 3': [35, 37, 7, 2, 19],
    'Question 4': [32, 11, 9, 15, 33],
    'Question 5': [21, 29, 5, 5, 40],
    'Question 6': [8, 19, 5, 30, 38]
}


def survey(results, category_names):
    """
    Parameters
    ----------
    results : dict
        A mapping from question labels to a list of answers per category.
        It is assumed all lists contain the same number of entries and that
        it matches the length of *category_names*.
    category_names : list of str
        The category labels.
    """
    labels = list(results.keys())
    data = np.array(list(results.values()))
    data_cum = data.cumsum(axis=1)
    category_colors = plt.get_cmap('RdGy_r')(
        np.linspace(0.15, 0.85, data.shape[1]))

    fig, ax = plt.subplots(figsize=(15, 7))
    ax.invert_yaxis()
    ax.xaxis.set_visible(False)
    ax.set_xlim(0, np.sum(data, axis=1).max())

    for i, (colname, color) in enumerate(zip(category_names, category_colors)):
        widths = data[:, i]
        starts = data_cum[:, i] - widths
        ax.barh(labels, widths, left=starts, height=0.5,
                label=colname, color=color)
        xcenters = starts + widths / 2

        r, g, b, _ = color
        text_color = 'white' if r * g * b < 0.5 else 'darkgrey'
        for y, (x, c) in enumerate(zip(xcenters, widths)):
            ax.text(x, y, str(int(c)), ha='center', va='center',
                    color=text_color)
    ax.legend(ncol=len(category_names), bbox_to_anchor=(0, 1),
              loc='lower left', fontsize='large')

    return fig, ax


survey(results, category_names)
Out[354]:
(<Figure size 1080x504 with 1 Axes>,
 <matplotlib.axes._subplots.AxesSubplot at 0x1af88838860>)
In [359]:
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(14)
y = np.sin(x / 2)

plt.figure(figsize=(15, 10))
plt.step(x, y + 2, label='pre (default)')
plt.plot(x, y + 2, 'C0o', alpha=0.5)

plt.step(x, y + 1, where='mid', label='mid')
plt.plot(x, y + 1, 'C1o', alpha=0.5)

plt.step(x, y, where='post', label='post')
plt.plot(x, y, 'C2o', alpha=0.5)

plt.legend(title='Parameter where:')
plt.show()
In [372]:
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 5))
plt.subplot(221, projection='aitoff')
plt.title('aitoff')

plt.subplot(222, projection='hammer')
plt.title('hammer')
plt.grid(True)

plt.subplot(223, projection='lambert')
plt.title('lambert')
plt.grid(True)

plt.subplot(224, projection='mollweide')
plt.title('mollweide')
plt.grid(True)
plt.show()
In [383]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 3
Y, X = np.mgrid[-w:w:100j, -w:w:100j]
U = -1 - X**2 + Y
V = 1 + X - Y**2
speed = np.sqrt(U**2 + V**2)

fig = plt.figure(figsize=(10, 13))
gs = gridspec.GridSpec(nrows=3, ncols=2, height_ratios=[1, 1, 2])

#  Varying density along a streamline
ax0 = fig.add_subplot(gs[0, 0])
ax0.streamplot(X, Y, U, V, density=[0.5, 1])
ax0.set_title('Varying Density')

ax1 = fig.add_subplot(gs[0, 1])
strm = ax1.streamplot(X, Y, U, V, color=U, linewidth=2, cmap='jet')
fig.colorbar(strm.lines)
ax1.set_title('Varying Color')

#  Varying line width along a streamline
ax2 = fig.add_subplot(gs[1, 0])
lw = 5*speed / speed.max()
ax2.streamplot(X, Y, U, V, density=0.6, color='k', linewidth=lw)
ax2.set_title('Varying Line Width')

# Controlling the starting points of the streamlines
seed_points = np.array([[-2, -1, 0, 1, 2, -1], [-2, -1,  0, 1, 2, 2]])

ax3 = fig.add_subplot(gs[1, 1])
strm = ax3.streamplot(X, Y, U, V, color=U, linewidth=2,
                     cmap='autumn', start_points=seed_points.T)
fig.colorbar(strm.lines)
ax3.set_title('Controlling Starting Points')

# Displaying the starting points with blue symbols.
ax3.plot(seed_points[0], seed_points[1], 'bo')
ax3.set(xlim=(-w, w), ylim=(-w, w))

# Create a mask
mask = np.zeros(U.shape, dtype=bool)
mask[40:60, 40:60] = True
U[:20, :20] = np.nan
U = np.ma.array(U, mask=mask)

ax4 = fig.add_subplot(gs[2, 0])
ax4.streamplot(X, Y, U, V, color='r')
ax4.set_title('Streamplot with Masking')

ax4.imshow(~mask, extent=(-w, w, -w, w), alpha=0.5,
          interpolation='nearest', cmap='gray', aspect='auto')
ax4.set_aspect('equal')

plt.tight_layout()
plt.show()
In [412]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

n = 20

Z = np.arange(n**2).reshape((n, n))
Z[:, 10:] = 1


fig=plt.figure()
plt.imshow(Z)
plt.imshow(Z)
plt.show()
In [392]:
import matplotlib.pyplot as plt
import numpy as np


def func3(x, y):
    return (1 - x / 2 + x**5 + y**3) * np.exp(-(x**2 + y**2))

# make these smaller to increase the resolution
dx, dy = 0.05, 0.05

x = np.arange(-3.0, 3.0, dx)
y = np.arange(-3.0, 3.0, dy)
X, Y = np.meshgrid(x, y)

extent = np.min(x), np.max(x), np.min(y), np.max(y)
fig = plt.figure(frameon=False, figsize=(10, 8))

Z1 = np.add.outer(range(8), range(8)) % 2  # chessboard
im1 = plt.imshow(Z1, cmap=plt.cm.gray, interpolation='nearest',
                 extent=extent)

Z2 = func3(X, Y)

im2 = plt.imshow(Z2, cmap=plt.cm.jet, alpha=.9, interpolation='bilinear',
                 extent=extent)

plt.show()
In [394]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.axes_rgb import RGBAxes


def get_demo_image():
    # prepare image
    delta = 0.5

    extent = (-3, 4, -4, 3)
    x = np.arange(-3.0, 4.001, delta)
    y = np.arange(-4.0, 3.001, delta)
    X, Y = np.meshgrid(x, y)
    Z1 = np.exp(-X**2 - Y**2)
    Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
    Z = (Z1 - Z2) * 2

    return Z, extent


def get_rgb():
    Z, extent = get_demo_image()

    Z[Z < 0] = 0.
    Z = Z / Z.max()

    R = Z[:13, :13]
    G = Z[2:, 2:]
    B = Z[:13, 2:]

    return R, G, B

fig = plt.figure(figsize=(10, 7))
ax = RGBAxes(fig, [0.1, 0.1, 0.8, 0.8])

r, g, b = get_rgb()
kwargs = dict(origin="lower", interpolation="nearest")
ax.imshow_rgb(r, g, b, **kwargs)

ax.RGB.set_xlim(0., 9.5)
ax.RGB.set_ylim(0.9, 10.6)

plt.show()
In [423]:
import numpy as np
import matplotlib.pyplot as plt


def koch_snowflake(order, scale=10):
    def _koch_snowflake_complex(order):
        if order == 0:
            # initial triangle
            angles = np.array([0, 120, 240]) + 90
            return scale / np.sqrt(3) * np.exp(np.deg2rad(angles) * 1j)
        else:
            ZR = 0.5 - 0.5j * np.sqrt(3) / 3

            p1 = _koch_snowflake_complex(order - 1)  # start points
            p2 = np.roll(p1, shift=-1)  # end points
            dp = p2 - p1  # connection vectors

            new_points = np.empty(len(p1) * 4, dtype=np.complex128)
            new_points[::4] = p1
            new_points[1::4] = p1 + dp / 3
            new_points[2::4] = p1 + dp * ZR
            new_points[3::4] = p1 + dp / 3 * 2
            return new_points

    points = _koch_snowflake_complex(order)
    x, y = points.real, points.imag
    return x, y

colors = ['r', 'g', 'b', 'k', 'y', 'c', 'm']
plt.figure(figsize=(15, 5))
for i in range(3):
    index1 = int(np.random.random()*7)
    index2 = int(np.random.random()*7)
    plt.subplot(1, 3, i+1)
    x, y = koch_snowflake(2+i)
    plt.fill(x, y, facecolor=colors[index1], edgecolor=colors[index2], linewidth=3)
plt.show()
In [426]:
import matplotlib.pyplot as plt


def handle_close(evt):
    print('Closed Figure!')

fig = plt.figure()
fig.canvas.mpl_connect('close_event', handle_close)

plt.text(0.35, 0.5, 'Close Me!', dict(size=30))
plt.show()
In [2]:
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)


def format_axes(fig):
    for i, ax in enumerate(fig.axes):
        ax.text(0.5, 0.5, "ax%d" % (i+1), va="center", ha="center")
        ax.tick_params(labelbottom=False, labelleft=False)

fig = plt.figure(constrained_layout=True, figsize=(12, 10))

gs = GridSpec(3, 3, figure=fig)
ax1 = fig.add_subplot(gs[0, :])
ax1.plot(t, s, linewidth=2)
# identical to ax1 = plt.subplot(gs.new_subplotspec((0, 0), colspan=3))
ax2 = fig.add_subplot(gs[1, :-1])
ax2.plot(t, s, linewidth=2)
ax3 = fig.add_subplot(gs[1:, -1])
ax3.plot(t, s, linewidth=2)
ax4 = fig.add_subplot(gs[-1, 0])
ax4.plot(t, s, linewidth=2)
ax5 = fig.add_subplot(gs[-1, -2])
ax5.plot(t, s, linewidth=2)

fig.suptitle("GridSpec")
format_axes(fig)

plt.show()
In [39]:
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)


def format_axes(fig):
    for i, ax in enumerate(fig.axes):
        ax.text(0.5, 0.5, "ax%d" % (i+1), va="center", ha="center")
        ax.tick_params(labelbottom=False, labelleft=False)

fig = plt.figure(constrained_layout=True, figsize=(12, 10))

gs = GridSpec(3, 3, figure=fig)
plt.xkcd()
ax1 = fig.add_subplot(gs[0, :])
ax1.plot(t, s, linewidth=2)
# identical to ax1 = plt.subplot(gs.new_subplotspec((0, 0), colspan=3))
ax2 = fig.add_subplot(gs[1, :-1])
ax2.plot(t, s, linewidth=2)
ax3 = fig.add_subplot(gs[1:, -1])
ax3.plot(t, s, linewidth=2)
ax4 = fig.add_subplot(gs[-1, 0])
ax4.plot(t, s, linewidth=2)
ax5 = fig.add_subplot(gs[-1, -2])
ax5.plot(t, s, linewidth=2)

fig.suptitle("GridSpec")
format_axes(fig)

plt.show()
In [38]:
import numpy as np
import matplotlib.pyplot as plt


mean = [0, 0]
cov = [[1, 1], [1, 2]]
x, y = np.random.multivariate_normal(mean, cov, 3000).T
v = np.sin(x**2+y**2)

#设置坐标轴和网格配置方式
fig = plt.figure(figsize=(12, 8))
grid = plt.GridSpec(4,4, hspace=0.2, wspace=0.2)
main_ax = fig.add_subplot(grid[:-1, 1:])
y_hist = fig.add_subplot(grid[:-1, 0], xticklabels=[], sharey = main_ax)
x_hist = fig.add_subplot(grid[-1, 1:], yticklabels=[], sharex = main_ax)

#主坐标轴画散点图
main_ax.scatter(x, y, s=200, c=v, marker=r'$\spadesuit$', alpha=0.2, cmap='jet')

#次坐标轴画频次直方图
x_hist.hist(x, 40, histtype='stepfilled', orientation='vertical', facecolor='b', edgecolor='k')
x_hist.invert_yaxis()
y_hist.hist(y, 40, histtype='stepfilled', orientation='horizontal', facecolor='r', edgecolor='k')
y_hist.invert_xaxis()

plt.show()
In [20]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)

plt.figure(figsize=(12, 8))
ax1 = plt.axes()
plt.plot(t, s, linewidth=2, color='r')
plt.text(5, 0, 'ax1', fontsize=25)

ax2 = plt.axes([0.65, 0.65, 0.2, 0.2])
plt.plot(t, s, linewidth=2, color='r')
plt.text(2, -0.5, 'ax2', fontsize=15)

ax3 = plt.axes([0.2, 0.2, 0.2, 0.2])
plt.plot(t, s, linewidth=2, color='r')
plt.text(2, -0.5, 'ax3', fontsize=15)

ax4 = plt.axes([0.65, 0.2, 0.2, 0.2])
plt.plot(t, s, linewidth=2, color='r')
plt.text(2, -0.5, 'ax4', fontsize=15)

ax5 = plt.axes([0.2, 0.65, 0.2, 0.2])
plt.plot(t, s, linewidth=2, color='r')
plt.text(2, -0.5, 'ax5', fontsize=15)

plt.show()
In [22]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)

fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_axes([0.1, 0.5, 0.8, 0.4], xticklabels=[])
ax1.plot(t, s, linewidth=2, color='r')

ax2 = fig.add_axes([0.1, 0.1, 0.8, 0.4])
ax2.plot(t, s, linewidth=2, color='b')

plt.show()
In [29]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)

fig, ax = plt.subplots(2, 3, sharex='col', sharey='row', figsize=(12, 8))
for i in range(2):
    for j in range(3):
        ax[i, j].plot(t, s, linewidth=2, color='r')
        ax[i, j].text(2.5, 0.5, str((i, j)),fontsize=18, ha='center')
plt.show()
In [4]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s = np.sin(t**2) * np.exp(-(t-t[500])**2)

plt.style.use('ggplot')

fig, ax = plt.subplots(2, 3, sharex='col', sharey='row', figsize=(12, 8))
for i in range(2):
    for j in range(3):
        ax[i, j].plot(t, s, linewidth=2, color='g')
        ax[i, j].text(2.5, 0.5, str((i, j)),fontsize=18, ha='center')
plt.suptitle('ggplot style', fontsize=25)
plt.show()
In [433]:
import matplotlib.pyplot as plt
np.random.seed(1)
x = np.arange(5)
y = np.random.randn(5)

fig, axes = plt.subplots(ncols=2, figsize=(12, 5))

vert_bars = axes[0].bar(x, y, color='b', align='center')
horiz_bars = axes[1].barh(x, y, color='r', align='center')

axes[0].axhline(0, color='gray', linewidth=2)
axes[1].axvline(0, color='gray', linewidth=2)
plt.show()
In [6]:
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(19680801)

mu1, sigma1 = 100, 15
mu2, sigma2 = 80, 15
x1 = mu1 + sigma1 * np.random.randn(10000)
x2 = mu2 + sigma2 * np.random.randn(10000)

# the histogram of the data
# 50:将数据分成50组
# facecolor:颜色;alpha:透明度
# density:是密度而不是具体数值
plt.figure(figsize=(10, 8))
n1, bins1, patches1 = plt.hist(x1, 50, density=True, facecolor='g', alpha=1)
n2, bins2, patches2 = plt.hist(x2, 50, density=True, facecolor='r', alpha=0.2)

# n:概率值;bins:具体数值;patches:直方图对象。

plt.xlabel('Smarts')
plt.ylabel('Probability')
plt.title('Histogram of IQ')

plt.text(110, .025, r'$\mu=100,\ \sigma=15$')
plt.text(50, .025, r'$\mu=80,\ \sigma=15$')

# 设置x,y轴的具体范围
plt.axis([40, 160, 0, 0.03])
plt.grid(True)
plt.show()
In [435]:
import matplotlib.pyplot as plt
import numpy as np

np.random.seed(19680801)

n_bins = 10
x = np.random.randn(1000, 3)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
ax0, ax1, ax2, ax3 = axes.flatten()

colors = ['red', 'tan', 'lime']
ax0.hist(x, n_bins, density=True, histtype='bar', color=colors, label=colors)
ax0.legend(prop={'size': 10})
ax0.set_title('bars with legend')

ax1.hist(x, n_bins, density=True, histtype='barstacked')
ax1.set_title('stacked bar')

ax2.hist(x,  histtype='barstacked', rwidth=0.9)

ax3.hist(x[:, 0], rwidth=0.9)
ax3.set_title('different sample sizes')

fig.tight_layout()
plt.show()
In [9]:
import numpy as np
import matplotlib.pyplot as plt

size = 5
a = np.random.random(size)
b = np.random.random(size)
c = np.random.random(size)

x = np.arange(size)

# 这里使用的是偏移
plt.figure(figsize=(12, 9))
plt.bar(x, a, width=0.5, label='a',fc='r')
plt.bar(x, b, bottom=a, width=0.5, label='b', fc='g')
plt.bar(x, c, bottom=a+b, width=0.5, label='c', fc='b')

plt.ylim(0, 2.5)
plt.legend()
plt.grid(True)
plt.title('Stacking histograms', fontsize=25)
plt.show()
In [39]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
s1 = np.sin(t**2) * np.exp(-(t-t[500])**2) - 0.5
s2 = np.cos(t**2) * np.exp(-(t-t[500])**2)

plt.figure(figsize=(15, 10))
plt.subplot(211)
plt.plot(t, s1, 'g', linewidth=2)
plt.plot(t, s2, 'b', linewidth=2)
plt.fill_between(t, s1, s2, color='r')

plt.subplot(212)
plt.fill_betweenx(t, t[100], t[200], color='g')
plt.fill_betweenx(t, t[500], t[700], color='r')
plt.fill_between(t, 5, 7, color='b', alpha=0.5)
plt.xlim(t[0], t[-1])

plt.show()
In [90]:
import numpy as np
import matplotlib.pyplot as plt

n = 11
m = 10
t1 = 5
t = np.linspace(0, t1, n)
d = np.random.random((m, n))

plt.figure(figsize=(15, 6))
for i in range(10):
    if i < 1:
        plt.fill_between(t, np.zeros(n), d[i], label=str(i+1))
    else:            
        plt.fill_between(t, d[i-1], d[i], label=str(i+1))
plt.legend()
plt.xlim(0, 5.5)
plt.ylim(0, 1)
plt.show()
In [38]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import PolyCollection
from matplotlib.colors import colorConverter

import scipy.io as sio

t = np.linspace(0,10,2000)

z = np.zeros([5,2000])

for i in range(5):
    z[i,:] = np.sin(2*np.pi*t*(i+1)*0.2)
    
x = [np.linspace(0,9,5)]*2000
y = [np.linspace(0,1999,2000)]*5
x = np.array(x).T
y = np.array(y)

fig = plt.figure(figsize=(6,2.5))

ax = fig.gca(projection='3d')
plt.gca().patch.set_facecolor('white')


ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

for i in range(5):
    ax.plot(x[i,:],y[i,:],z[i,:])

ax.view_init(20,45)     
fig.set_size_inches(12, 7)
plt.show()
In [12]:
import matplotlib.pyplot as plt
import numpy as np

n = 201
t1 = 2

t = np.linspace(0, t1, n)
s = t**3 + np.random.random(n)
n1 = np.random.random(n) + s - 1.5
n2 = np.random.random(n) + s + 0.5

plt.figure(figsize=(15, 6))
plt.fill_between(t, n1, n2, facecolor='r', alpha=0.25)
plt.plot(t, s, 'r', linewidth=2)
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()
In [68]:
import numpy as np
import matplotlib.pyplot as plt

nl = 301
t1 = 3
ns = 250
t = np.linspace(0, 3, nl)
s = np.array([np.sin(2*np.pi*(t[i])) if i >=  150 and i < 200 else 0 for i in range(nl)])
#s = np.sin(np.pi*t)
t = t[: ns]
n = 75

plt.figure(figsize=(6, 8))
for i in range(n):
    index = int(np.random.random() * 40)
    ss = abs(s[index: index+ns]) + i
    ss2 = np.ones(ns)*i
    #plt.plot(t, ss+i)
    #plt.plot(t, ss2, 'k', linewidth=2)
    plt.fill_between(t, ss, ss2, facecolor='b', edgecolor='k', linewidth=1)
    
plt.show()
In [15]:
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 1001)
n1 = np.random.random(1001)
n2 = np.random.random(1001)
e = np.exp(-t)
en1 = e + (n1-n1.mean())/5.
en2 = e + (n2-n2.mean())/5.
e1 = en1 - en1.max()/5
e2 = en1 + en1.max()/5
plt.figure(figsize=(15,10))
plt.fill_between(t, e1, e2, facecolor='green', alpha=0.5)
plt.subplot(111)
plt.plot(t, en1, 'b', linewidth=1, label=r'$Noised1 : e^{-t}$')
plt.subplot(111)
plt.plot(t, en2, 'r', linewidth=1, label=r'$Noised2 : e^{-t} $')
plt.legend(loc='best', fontsize=20)
plt.xlim(0, 10)
plt.grid()
plt.show()
D:\Anaconda\lib\site-packages\matplotlib\cbook\deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [16]:
import numpy as np
import matplotlib.pyplot as plt

# 设置每环的宽度
size = 0.3
vals = np.array([[60., 32.], [37., 40.], [29., 10.]])

# 通过get_cmap随机获取颜色
cmap = plt.get_cmap("tab20c")
outer_colors = cmap(np.arange(3)*4)
inner_colors = cmap(np.array([1, 2, 5, 6, 9, 10]))

print(vals.sum(axis=1))
# [92. 77. 39.]

plt.figure(figsize=(12, 9))
plt.pie(vals.sum(axis=1), radius=1, colors=outer_colors,
       wedgeprops=dict(width=size, edgecolor='w'))
print(vals.flatten())
# [60. 32. 37. 40. 29. 10.]

plt.pie(vals.flatten(), radius=1-size, colors=inner_colors,
       wedgeprops=dict(width=size, edgecolor='w'))

# equal 使得为正圆
plt.axis('equal') 
plt.show()
[92. 77. 39.]
[60. 32. 37. 40. 29. 10.]
In [16]:
import numpy as np
import matplotlib.pyplot as plt

labels = 'Frogs', 'Hogs', 'Dogs', 'Logs'
sizes = [15, 30, 45, 10]
explode = (0, 0.1, 0, 0)  # only "explode" the 2nd slice (i.e. 'Hogs')

fig1, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 8))
ax1.pie(sizes, labels=labels, autopct='%1.1f%%', shadow=True)
ax1.axis('equal')
ax2.pie(sizes, autopct='%1.2f%%', shadow=True, startangle=90, explode=explode,
    pctdistance=1.12)
ax2.axis('equal')
ax2.legend(labels=labels, loc='upper right')

plt.show()
In [19]:
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(19680801)

N = 10
theta = np.linspace(0.0, 2 * np.pi, N, endpoint=False)
radii = 10 * np.random.rand(N)
width = np.pi / 4 * np.random.rand(N)

plt.figure(figsize=(12, 9))
ax = plt.subplot(111, projection='polar')
bars = ax.bar(theta, radii, width=width, bottom=0.0)
# left表示从哪开始,
# radii表示从中心点向边缘绘制的长度(半径)
# width表示末端的弧长

# 自定义颜色和不透明度
for r, bar in zip(radii, bars):
    bar.set_facecolor(plt.cm.jet(r / 10.))
    bar.set_alpha(0.9)

plt.show()
In [7]:
import numpy as np
import matplotlib.pyplot as plt


data = np.random.rand(10)
data2 = np.random.rand(10, 10)
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 7))
ax1.boxplot(data)
ax2.boxplot(data2, vert=False)

plt.show()
In [28]:
import matplotlib.pyplot as plt
import matplotlib as mpl

fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)

cmap = mpl.cm.cool
norm = mpl.colors.Normalize(vmin=5, vmax=10)

cb1 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
                                norm=norm,
                                orientation='horizontal')
cb1.set_label('Some Units')
fig.show()
D:\Softwares\Anaconda\lib\site-packages\matplotlib\figure.py:457: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
In [37]:
# import matplotlib.pyplot as plt
# import numpy as np
# from matplotlib import animation
 
# fig, ax = plt.subplots()
 
# x = np.arange(0, 4 * np.pi, 0.01)
# line, = ax.plot(x, np.sin(x))
 
 
# def animate(i):
#     line.set_ydata(np.sin(x + i / 200))
#     return line,
 
 
# def init():
#     line.set_ydata(np.sin(x))
#     return line,
 
 
# ani = animation.FuncAnimation(fig=fig, func=animate, frames=100,
#                               init_func=init, interval=20, blit=False)
# plt.show()
In [8]:
import matplotlib.pyplot as plt
import matplotlib as mpl

cmap = mpl.colors.Colormap('jet')
print(cmap.N)
256
In [15]:
from obspy.clients.fdsn import Client
from obspy import UTCDateTime

client = Client('IRIS')
lt = 2500
starttime = UTCDateTime('2017-12-15T16:47:56.000')
st = client.get_waveforms('IC', 'XAN', '*', 'BH?', starttime, starttime+lt)
st.filter('bandpass', freqmin=0.005, freqmax=0.05, zerophase=True)
st.plot(color='r', linewidth=0.75, alpha=0.2)
Out[15]:
In [20]:
from obspy.clients.fdsn import Client
from obspy import UTCDateTime

client = Client('IRIS')
lt = 1500

t = UTCDateTime("2017-11-12T18:18:19.000")
st = client.get_waveforms('IC', 'XAN', '*', 'BH?', t+300, t+lt)
st.filter('bandpass', freqmin=0.005, freqmax=0.1, zerophase=True)
st.plot(color='r', linewidth=0.75, alpha=0.2)
Out[20]:
In [2]:
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
import numpy as np
import matplotlib.pyplot as plt

starttime = UTCDateTime('2019-08-01T00:00:00.000')
endtime = UTCDateTime('2019-08-20T23:59:59.999')
minlongitude=100
maxlongitude=140
minlatitude=10
maxlatitude=45

client = Client('IRIS')
inventory = client.get_stations(minlongitude=minlongitude,
                               maxlongitude=maxlongitude,
                               minlatitude=minlatitude,
                               maxlatitude=maxlatitude,
                               starttime=starttime,
                               endtime=endtime)
#print(inventory)
coord = list()
with open('sta.txt', 'w') as fin:
    for net in inventory:
        for sta in net:
            coord.append([sta.longitude, sta.latitude, sta.elevation])
            out = ' '.join([net.code, sta.code, str(sta.longitude), str(sta.latitude), str(sta.elevation), '\n'])
            fin.write(out)
    #         print(net.code)
    #         print(sta.code)
    #         print(sta.latitude)
    #         print(sta.longitude)
    #         print(sta.elevation)
    #         print('-----------------------')
     
coord = np.array(coord)
center = np.mean(coord, axis=0)
print(center)
dist = ((coord[:, 0]-center[0])**2+(coord[:, 1]-center[1])**2)**0.5
plt.figure(figsize=(9, 8.5))
plt.scatter(coord[:, 0], coord[:, 1], c=dist, s=200, marker='v', cmap='jet')
plt.xlim(minlongitude-maxlongitude*0.1, maxlongitude+maxlongitude*0.1)
plt.ylim(minlatitude-maxlatitude*0.1, maxlatitude+maxlatitude*0.1)
plt.xlabel('Longitude [$^o$]', fontsize=20)
plt.ylabel('Latitude [$^o$]', fontsize=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.title('Get station info. from IRIS', fontsize=25)
plt.grid(True)
cbar=plt.colorbar(extend='both')
cbar.set_label('Distance [$^o$]')
plt.show()
print(len(coord))
[ 118.173702     32.8960079   546.37314488]
283
In [81]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111, projection='3d')

X = np.random.random(4)
Y = np.random.random(4)
Z = np.random.random(4)

ax.plot_trisurf(X, Y, Z)
plt.show()
In [80]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

data = np.random.randint(0, 255, size=[40, 40, 40])

x, y, z = data[0], data[1], data[2]
plt.figure(figsize=(7, 7))
ax = plt.subplot(projection='3d')  # 创建一个三维的绘图工程
#  将数据点分成三部分画,在颜色上有区分度
ax.scatter(x[:10], y[:10], z[:10], c='y', s=100)  # 绘制数据点
ax.scatter(x[10:20], y[10:20], z[10:20], c='r', s=100)
ax.scatter(x[30:40], y[30:40], z[30:40], c='g', s=100)

ax.set_zlabel('Z')  # 坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')

fig=plt.figure(figsize=(7, 7))
ax = Axes3D(fig)
X = np.arange(-4, 4, 0.25)
Y = np.arange(-4, 4, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)

# 具体函数方法可用 help(function) 查看,如:help(ax.plot_surface)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
plt.show()
In [75]:
import numpy as np

with open('sta.txt') as fin:
    print(fin.readlines()[0])
AD DLV 108.481499 11.952 50.0 

In [50]:
from pylab import *
eqs = []
eqs.append((r'$\int_{-\infty}^{\infty}f(x)e^{-i\omega x}dx $'))
axes([0.025, 0.025, 0.95, 0.95])
text(0.5, 0.85, eqs[0], ha='center', va='center', color='r', alpha=1, transform=gca().transAxes,
fontsize=20, clip_on=True)
text(0.5, 0.6, eqs[0], ha='center', va='center', color='b', alpha=1, transform=gca().transAxes,
fontsize=20, clip_on=True)
text(0.5, 0.4, eqs[0], ha='center', va='center', color='k', alpha=1, transform=gca().transAxes,
fontsize=20, clip_on=True)
text(0.5, 0.15, r'$F(k) = \sum_{n=0}^{N-1}x(n)e^{\frac{-i \pi 2nk}{N}}$', ha='center', va='center',\
color='k', alpha=1, transform=gca().transAxes, fontsize=25, clip_on=True)
xticks([])
yticks([])
title('Mathematical equations', fontsize=25)
show()
#savefig('eqs.pdf
In [4]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

plt.figure(figsize=(20, 5.5))
plt.subplot(131)
map = Basemap(projection='cyl', lat_0=0, lon_0=0)
map.drawcoastlines(color='blue')
plt.title('coastlines', fontsize=20)
plt.subplot(132)
map = Basemap()
map.drawcoastlines(color='black')
map.drawcountries(color='red')
plt.title('coastlines and countries\' lines', fontsize=20)
plt.subplot(133)
map=Basemap()
map.drawcoastlines(color='black')
map.drawstates(color='green')
map.drawrivers(color='blue')
plt.title('coastlines, states\' lines and rivers', fontsize=20)
plt.suptitle('Simple geoMap', fontsize=35)
plt.show()
In [6]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
map = Basemap(projection = 'cyl')
map.drawmapboundary(fill_color = 'aqua')
map.fillcontinents(color = 'coral', lake_color = 'aqua')
map.drawcoastlines()
plt.title('fill continents and oceans', fontsize=30)
plt.show()
In [8]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

# setup Lambert Conformal basemap.
# set resolution=None to skip processing of boundary datasets.
plt.figure(figsize=(20, 6))
plt.subplot(131)
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution=None,lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.bluemarble(scale=0.5)
plt.title('m.bluemarble()', fontsize=20)
plt.subplot(132)
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution=None,lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.shadedrelief(scale=0.5)
plt.title('m.shadedrelief()', fontsize=20)
plt.subplot(133)
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution=None,lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.etopo(scale=0.5)
plt.title('m.etopo()', fontsize=20)
plt.suptitle('Basemap', fontsize=30)
plt.xticks(fontsize=20)
plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
In [10]:
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 8))
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution=None,lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.bluemarble()
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=20)
meridians = np.arange(180.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=20)
plt.title('setting ticks of geoMap', fontsize=30)
plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
In [12]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
map = Basemap(llcrnrlon=3.75,llcrnrlat=39.75,urcrnrlon=4.35,urcrnrlat=40.15, epsg=5520)
map.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 1500, verbose= True)
plt.title('Get image from ArcGIS', fontsize=30)
plt.show()
http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_Imagery_World_2D/MapServer/export?bbox=1564270.9620172689,4401598.726055895,1615017.0560379555,4446612.184939481&bboxSR=5520&imageSR=5520&size=1500,1330&dpi=96&format=png32&transparent=true&f=image
In [38]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(20, 13))
plt.subplot(341)
map = Basemap(projection='ortho', lat_0=0, lon_0=0)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral', lake_color='aqua')
map.drawcoastlines(color='black')
plt.title('ortho projection', fontsize=20)

plt.subplot(342)
map = Basemap(projection='hammer', lat_0=0, lon_0=0)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral', lake_color='aqua')
map.drawcoastlines(color='black')
plt.title('hammer projection', fontsize=20)

plt.subplot(343)
map = Basemap(projection = 'aeqd', lon_0 = 10, lat_0 = 50)
map.drawmapboundary(fill_color = 'aqua')
map.fillcontinents(color = 'coral', lake_color = 'aqua')
map.drawcoastlines()
plt.title('aeqd projection', fontsize=20)

plt.subplot(344)
m = Basemap(width=28000000, height=28000000, projection='gnom', lat_0=0, lon_0=0)
m.drawmapboundary(fill_color='white')
m.drawcoastlines(linewidth=0.5)
m.fillcontinents(color='gray',lake_color='white')
m.drawparallels(np.arange(-80,81,20))
m.drawmeridians(np.arange(-180,180,20))
plt.title('gnom Projection', fontsize=20)

plt.subplot(345)
m = Basemap(projection='moll',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,420.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title('moll Projection', fontsize=20)

plt.subplot(346)
m = Basemap(projection='robin',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("robin Projection", fontsize=20)

plt.subplot(347)
m = Basemap(projection='eck4',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("eck4 Projection", fontsize=20)

plt.subplot(348)
m = Basemap(projection='kav7',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("kav7 Projection", fontsize=20)

plt.subplot(349)
m = Basemap(projection='mbtfpq',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("mbtfpq Projection", fontsize=20)

plt.subplot(3, 4, 10)
m = Basemap(projection='sinu',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("sinu projection", fontsize=20)

plt.subplot(3, 4, 11)
m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80,\
llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,91.,30.))
m.drawmeridians(np.arange(-180.,181.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("merc Projection", fontsize=20)
plt.show()
In [16]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(9, 6))
map = Basemap(projection='ortho', lat_0=32, lon_0=120)
map.bluemarble()
plt.title('ortho projection and \'bluemarble\' cpt', fontsize=30)
plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
In [18]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(5, 6))
map = Basemap(llcrnrlon = 119.3, llcrnrlat = 20.7, urcrnrlon = 124.6, urcrnrlat = 26,
resolution = 'h', epsg = 3415)
map.drawmapboundary(fill_color = 'aqua')
map.fillcontinents(color = '#003300', lake_color = '#0000FF')
map.drawcoastlines()
plt.show()
In [20]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
sta_lat = np.linspace(30,50, 5)
sta_lon = np.linspace(-120,-80,5)
plt.figure(figsize=(12, 8))
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution=None,lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.shadedrelief()
for i, j in zip(sta_lon, sta_lat):
    x, y = m(i, j)
    m.scatter(x, y, marker='v', s=400, color='k')
    x, y = m(-70, 10)
    m.scatter(x, y, marker='*', s=400, color='r')
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=20)
meridians = np.arange(180.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=20)
plt.title('epicenter and stations', fontsize=30)
plt.show()
In [21]:
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from datetime import datetime
plt.figure(figsize=(10, 6))
# miller projection
map = Basemap(projection='mill',lon_0=180)
# plot coastlines, draw label meridians and parallels.
map.drawcoastlines()
map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0], fontsize=20)
map.drawmeridians(np.arange(map.lonmin,map.lonmax+30,60),labels=[0,0,0,1], fontsize=20)
# fill continents 'coral' (with zorder=0), color wet areas 'aqua'
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
# shade the night areas, with alpha transparency so the
# map shows through. Use current time in UTC.
date = datetime.utcnow()
CS=map.nightshade(date)
plt.title('Day/Night Map for %s (UTC)' % date.strftime("%d %b %Y %H:%M:%S"), fontsize=30)
plt.show()
In [43]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
lon = np.linspace(-100, 100, 21)
lat = np.zeros(21)
plt.figure(figsize=(9, 6))
map = Basemap(projection='hammer', lat_0=0, lon_0=0)
map.drawmapboundary(fill_color='aqua')
#map.fillcontinents(color='gray', lake_color='aqua')
map.drawcoastlines(color='black')
x, y = map(lon, lat)
map.scatter(x, y, marker='o', s=200, color='r')
x, y = map(-100, 0)
plt.text(x-5000000, y+500000, 'station', fontsize=25, color='blue')
plt.title('plt.text()', fontsize=25)
plt.show()
In [48]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
lon = np.random.random(21)*100-50
lat = np.random.random(21)*80-40
x = -100
y = 30
plt.figure(figsize=(12, 9))
map = Basemap(projection='kav7', lon_0=0, resolution='c')
map.drawcoastlines(color='k', linewidth=1)
#map.fillcontinents(color='#EEEEEE', lake_color='aqua')
for i, j in zip(lon, lat):
    map.drawgreatcircle(i, j, x, y, del_s=111.195, linewidth=2)
i, j = map(lon, lat)
m, n = map(x, y)
map.scatter(i, j, marker='v', s=300, color='k')
map.scatter(m, n, marker='*', s=500, color='red')
'''
x = np.linspace(-50, 50, 11)
y = np.linspace(-30, 30, 21)
X, Y = np.meshgrid(x, y)
v = np.random.random((21,11))
plt.figure(figsize=(8,6))
map.pcolormesh(X, Y, v, cmap='hot')
map.colorbar(extend='both')
'''
parallels = np.linspace(-90, 90, 7)
meridians = np.linspace(0, 360, 7)
map.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=16)
map.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=16)
plt.title('plot great circle arc', fontsize=25)
plt.show()
In [52]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# set up orthographic map projection with
# perspective of satellite looking down at 50N, 100W.
# use low resolution coastlines.
plt.figure(figsize=(20, 7.5))
plt.subplot(121)
map = Basemap(projection='ortho',lat_0=45,lon_0=-100,resolution='l')
# draw coastlines, country boundaries, fill continents.
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color='#EEEEEE',lake_color='aqua')
# draw the edge of the map projection region (the projection limb)
map.drawmapboundary(fill_color='aqua')
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
# make up some data on a regular lat/lon grid.
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# compute native map projection coordinates of lat/lon grid.
x, y = map(lons*180./np.pi, lats*180./np.pi)
# contour data over the map.
map.contourf(x,y,wave+mean, cmap='viridis')
plt.colorbar(extend='both', shrink=0.75)

plt.subplot(122)
map = Basemap(projection='ortho',lat_0=45,lon_0=-100,resolution='l')
# draw coastlines, country boundaries, fill continents.
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color='#EEEEEE',lake_color='aqua')
# draw the edge of the map projection region (the projection limb)
map.drawmapboundary(fill_color='aqua')
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
# make up some data on a regular lat/lon grid.
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# compute native map projection coordinates of lat/lon grid.
x, y = map(lons*180./np.pi, lats*180./np.pi)
# contour data over the map.
map.contour(x,y,wave+mean, cmap='jet')
plt.suptitle('contour lines over filled continent background', fontsize=20)
plt.show()
In [55]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
fig=plt.figure(figsize=(13, 6))
map = Basemap(projection='robin',
lat_0=0, lon_0=0)
lons = np.arange(-180, 190, 30)
lats = np.arange(-90, 100, 15)
data = np.indices((lats.shape[0], lons.shape[0]))
data = data[0] + data[1]
llons, llats = np.meshgrid(lons, lats)
ax = fig.add_subplot(121)
ax.set_title('Without transform_scalar', fontsize=20)
x, y = map(llons, llats)
map.pcolormesh(x, y, data, cmap='jet')
map.drawcoastlines()
ax = fig.add_subplot(122)
ax.set_title('Applying transform_scalar', fontsize=20)
data_interp, x, y = map.transform_scalar(data, lons, lats, 50, 30, returnxy=True, masked=True)
map.pcolormesh(x, y, data_interp, cmap='jet')
map.drawcoastlines()
plt.show()
In [57]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
x = np.linspace(-180, 180, 361)
y = np.linspace(-90, 90, 181)
X, Y = np.meshgrid(x, y)
v = np.sqrt(X**2+Y**2)
plt.figure(figsize=(10, 7))
map = Basemap(projection='hammer', lat_0=0, lon_0=0)
x, y = map(X, Y)
map.pcolormesh(x, y, v, cmap='RdBu_r')
map.drawcoastlines(color='k', linewidth=1)
map.drawmapboundary(color='k')
plt.colorbar(extend='both', shrink=0.45)
plt.title('Extra data plotted on basic map', fontsize=25)
plt.show()
In [59]:
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
x = np.linspace(-90, 90, 91)
y = np.linspace(-90, 90, 91)
X, Y = np.meshgrid(x,y)
v = np.exp(-X**2/10-Y**2/2) + np.sqrt(X**2+Y**2) #+ np.random.random((101, 121))*10
plt.figure(figsize=(10, 7))
map = Basemap(projection='ortho', lat_0=0, lon_0=0)
x, y = map(X, Y)
map.pcolormesh(x, y, v, cmap='nipy_spectral')
map.drawcoastlines(linewidth=1, color='k')
map.fillcontinents(color='#CCCCCC')
plt.colorbar(extend='both', shrink=0.5)
plt.title('Plot extra data on basemap', fontsize=20)
plt.show()